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ABSTRACT 
A  pair  or  discrete-time  queues  is  coupled  by  the  fact  that  the  arrivals  on  the 
trwo  are  not  independent.  We  derive  a  runctional  equational  for  the  bivariate 
generating  function  of  the  joint  distribution  of  queue  lengths  at  equilibrium. 
This  is  solved,  and  the  result  is  used  to  And  exact  expressions  for  the  mean 
and  variance  of  the  sum  of  the  queue  lengths  and  to  covariance  of  the  pair  of 
the  queue  lengths.  We  also  find  asymptotic  formulae  for  the  equilibrium  pro- 
babilities pij  as  ij-*oo  ,  as  well  as  some  asymtotic  results  when  the  arrivval 
rate  approaches  0  or  1. 


Introduction. 

The  following  discrete-time  queueing  problem  arises  in  the  design  of  parallel  computer 
networks  (see  [5]  for  details).  There  are  two  queues  and  two  independent  input  lines.  A 
customer  arriving  on  either  line  can  choose  to  join  either  queue  with  probability  1/2.  The 
arrivals  and  service  occur  at  discrete  time  steps.  At  each  step  there  is  a  probability  p  of  an 
arrival  on  each  line,  and  one  customer  on  each  non-empty  queue  is  served  (fig.  1).  Therefore 
at  each  step  0,  1  or  2  customers  can  arrive  on  a  given  queue,  but  the  arrivals  on  the  two 
queues  are  not  independent. 
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More  concretely,  there  are  two  types  of  customers,  A  and  B,  each  serviced  by  its  own 
queue,  each  input  line  carries  both  types,  and  each  customer  is  equally  likely  to  be  of  either 
type.  In  the  network  described  in  [5],  the  customers  are  memory  requests  and  A  and  B 
label  two  segments  of  memory. 

Letting  (I,  J)  denote  the  lengths  of  the  two  queues,  this  process  can  be  viewed  as  a 
random  walk  on  the  non-negative  integer  lattice  with  reflecting  boundaries,  whose  transition 
probabilities  are  shown  in  figure  2. 


This  is  an  irreducible,  aperiodic  Markov  chain.   It  will  be  shown  below  that  it  is  ergodic. 
There  is  therefore  a  unique  equilibrium  distribution 
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P{I  =1,   J  =  j) 


To  find  this  distribution  we  define  the  generating  function 
(0.1a) 
(0.1b) 


i,j  >0 

which  is  analytic  for  |x|,  \y\  <  1  and  continuous  for  |i|,  \y\  <  1 


The  equilibrium  equations  for  p,j  lead  to  a  functional  equation  for  p{x,  y),  whose  solution, 
subject  to  conditions  (0.1b),  provides  a  nice  illustration  of  two  distinct  apijrDUches  ])rcviousIy 
used  to  study  paired  continuous-time  queues:  the  boundary-value-problem  formulation  of 
Boxma  and  Cohen  [3],  and  the  automorphy  method  of  Flatto  and  McKean  [4]. 

In  Section  1  we  derive  the  functional  equation.  In  Section  2  we  derive  and  solve  the 
associated  boundary  value  problem,  obtaining  an  explicit  solution  in  terms  of  Jacobi  elliptic 


functions.  In  Section  3  we  show  how  the  automorphy  method  produces  independently  a 
solution  for  p{x,y)  in  terms  of  a  partial-fraction  expansion.  Section  4  presents  exact  results 
for  the  mean  and  variance  of  the  total  queue  size,  asymptotic  results  on  pj_,  for  i,  j  — *  oo, 
and  Limiting  behavior  for  the  arrival  rate  p  — >  0  and  p  — >  1. 

Percus  and  Percus   [7]  have  given  a  much  fuUer   analysis  of  the  marginal  probabilities 
associated  with  a  single  one  of  the  two  queues,  including  properties  of  the  output  process. 

1.   The  functional  equation. 

Let  T  be  the  transition  matrix  for  the  Markov  chain: 

The  equilibrium  probabilities  p  =  (Pij)  satisfy  T'p  =  p,  where  T*  is  the  adjoint  of  T.  For 
a  =  {a,j)  e  £i,  b  =  (6,j)  e  ^oo,  let  (a,  b)  =  J2,,j>o°-v^^j  ■  Then  p(x,y)  =  (p,z)  where 
z,j  :=  x't/'  ,  and  T'p  =  p  implies 

(1.1)  (p,z)  =  (r'p,z)  =  (p,Tz) 

Now 

(Tz).,  =   X^  P({t,j)^[k,l))x'y' 

k,l>0 


=   K3P((j,j)-.(z  +  n,fc  +  m))x"y-jxV 

\n,m  / 


where  the  homogeneity  of  the  walk  implies  that  R,j(x,y)  =  R(x,y)  is  independent  of  i,j 
for  i,  j  >  1.  Therefore,  treating  separately  the  t  =  0  and  j  =  0  terms,  (1.1)  gives 

p(i,  y)  =  R(x,  y)p{x,  y)  +  boundary  terms 

and  using  the  transition  probabilities  (fig.    2),  a  tedious  but  simple  calculation  gives  the 
functional  equation  for  p(x,y): 

(1.2)  (xy  -  r(x,  y))p{x,  j/)  =  (y  -  l)7-(x,  0)p(x,  0)  +  (x  -  l)7-(0,  y)p(0,  y) 

+  (x-l)(y-l)r(0,0)p(0,0) 

where 

•> 
(1.2a)  r{x,y)  =  xyR{x,y)  =  (l  -p+  |(^  +  y))" 

Note  that  I  —  (ij/)~^  r(i,  y)  is  simply  the  generating  function  of  the  homogeneous  dif- 
ference operator  satisfied  by  p  in  the  interior,  and  the  right  hand  side  of  (1.2)  represents 
boundary  terms.  The  remarkably  simple  form  of  the  boundary  terms  is  no  accident;  in  fact 
any  homogeneous  nearest-neighbor  random  walk  in  the  quarter-plane  wLU  satify  (12)  with 
appropriate  r(i,y). 

To  see  this,  associate  to  the  transition  diagram  at  point  {i,  j)  the  generating  function  of 
the  corresponding  difference  operator 

1 


and  let  r,,(x,y)  =  iyi?.;(x,y).  For  example,  the  diagram 


i- Vx 

has  r(i,y)  =  ;a;y(i-^  +  X +  y"^  +y)- 

For  a  homogeneous  walk,  r,(x,y)  is  mdependent  of.,,  for  z,;  >  1,  and  its  diagram 


has  r-(x,  y)  =  y;k,J=o  <^*'^    «  ■  .  ■  _  n    '  -^  i 

The  boundary  terms  are,  in  obvious  notation,  on  t  -  U,  j  >  i. 


corresponding  to 

rv{x,  y)  =  ^(x,  y)  -  r(0,  y)  +  xr(0,  y)  =  r(x,  y)  +  (x  -  l)r(0,  y) 

Similarly  we  get  on  j  =  0,  i  >  1: 

and  composing  these  operations  gives  for  1  =  j  =  0:  ^^ 

'•o(x,y)  =  {Tv)H(x,y) 

=  r(x,  y)  +  (x  -  l)r(0,  y)  +  (y  -  l)'-(0,  x) 
+  (x-l)(y-l)r(0,0) 

Therefore  we  get 

iyp(x,y)  =  y^  p.j'-.j(x,y)xV 

=  r(x,y)  y  p.,x'y^  +(rv{x,y)~r{x,y])Y.Po,y'  +{rH{x,y)-r{x,y))J2p.ox' 
+  (T-o(x,y)  -rv{x,y)  -rH(x,y)  +r(x,y))poo 

which  gives  equation  (1.2)  

Note  that  if  an  equilibrium  exists  then  setting  x  or  y  =  1  gives  the  generating  function 
for  the  marginal  distribution  of  the  length  of  a  single  queue: 


(1.3) 


p(x,  1)  =  Ci 
p(l,I/)  =  ^: 


X  -  1 


X  —  r(i, 1) 

y-1 


'y-  ''(i,y) 

where  the  normalization  p(l,  1)  =  1  gives  Ci  =  1  -  r^(l,  1),  C2  =  1  -  ry(l,  1). 


Now  r(i,  1)  =  Ao  +  AiX  +  A2X-  where  A,  =  I^J=o  °'J  •  Since  1  is  a  root  of  i  -r(x,  1),  the 
other  root  is  A0/A2;  similarly  the  roots  of  a;  — r(l,  y)  are  1  and  A0/A2,  where  A,  =  X]j=o°J'- 
Therefore  in  order  for  an  equilibrium  distribution  to  exist  it  is  necessary  that 

(1.4)  Ao  >  A2     and     Aq  >  A2 

These  are  just  the  well-known  conditions  for  ergodicity  of  the  associated  one-dimensional 
walks,  obtained  by  projecting  onto  I  ov  J  respectively.  Note  that  these  projections  are  them- 
selves Markov  chains,  since  J2m  ^((^'j)  ~^  C'"^))  i^  independent  of  j,  and  V],P((i,  j)  — » 
(I,  m))  is  independent  of  i.  But  this  means  that  (1.4)  is  also  sufficient  for  the  (7,  J)  chain  to 
be  ergodic,  since  if  each  projection  is  ergodic  then  for  fixed  Jq,  Jo,  and  A'^  sufficiently  large 
we  have 

Um   P(In  <  N  \Io)>  I  and    Um    P(J„  <  iV  |  Jo)  >  - 

where  (/„,  Jn)  is  the  position  after  n  steps.   Therefore 

hm   P(J„  <  TV,   Jn<N  \  (7o,  Jo))  >  0 

n— •  00 

which  implies  ergodicity  [2]. 

In  particular,  the  chain  defined  by  fig.  2  is  ergodic  for  p  <  1,  since  ,4o  =  ^0  =  (1  —  p/2)^ 
and  A2  =  A2  =  (p/2)'. 

2.   Solution  of  the  functional  equation:   Conformal  mapping. 

Note  that  setting  i  or  y  =  0  in  (1.2)  yields  an  identity,  so  (1.2)  alone  imposes  no  re- 
strictions on  p(x,0)  and  p(0,  y).  However  (0.1b)  implies  that  the  RHS  of  (1.2)  must  vanish 
whenever  xy  —  r(i,  y)  =  0  with  |x|,  \y\  <  1,  and  this  turns  out  to  be  sufficient  to  determine 
p(i,  0)  and  p(0,  y),  and  hence  p(i,  y),  uniquely  up  to  a  multiphcative  constant.  The  constant 
comes  from  the  normalization  p(l,  1)  =  1. 

Let  S  denote  the  complex  curve  xy  —  r{x,y)  =  0,  and  D^  —  {|x|  <  1,  \y\  <  l}.  Define 

.t(x)=^:i^^^^  +  1.(0,  0)p(0,0)  and 
a:  —  1  2 

,   ,        7-(0,y)p(0,y)        1 

g2(y)  =  :, +  -r(0,0)p(0,0) 

y-  1  2 

Dividing  through  (1.2)  by  [x  -  l)(y  -  1)  we  see  that  the  condition  that  the  RHS  of  (1.2) 
vanish  on  ^  n  D"  is  equivalent  to  finding  functions  g^  such  that 

(2.1)  9i(i)  +  92(y)  =0  for  {x,y)eSnD-,  x,y^  1 
with 

g,{z)  analytic  for  \z\  <  1,  continuous  for  |.:|  <  1,   -  ^  1,  and 
(2.1a)  ii5i(^  -  l)ffi(^)  =  '"(^'  0)P(1'  0)  =  1  -  ry(l,  1) 

Um(y  -  l)3o(y)  =  r(0,  l)p(0,  1)  =  1  -  r,(l,  1) 

y— 1 

where  p(l,  0)  and  p(0,  1)  have  been  obtained  from  (1-3).  Using  (1.2),  p(x,  y)  is  now  recovered 

by 

(2.2)  P(x,y)=  ^^^-^^^^^(gi(x)  +  g,(y)) 

xy  -  r(x,y) 
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From  now  on,  we  confine  ourselves  to  the  random  walk  given  by  fig.  2.  This  leads  to  two 
important  simplifications:  First,  p(x,y)  is  now  synimetric,  so  pi  =  g^-  Second,  r[x,y)  is  a 
quadratic  form,  which  makes  S  particularly  simple.  The  method  of  Cohen  and  Boxma[3] 
will  also  solve  the  more  general  case  (without  symmetry  and  with  more  general  r),  but  in 
much  less  explicit  form.  Flatto  and  McKean  solve  a  problem  with  a  similar  r,  but  without 
the  requirement  of  symmetry. 

Since  r(x,y)  (which  is  given  by  (1.2a))  and  p(x,  j/)  are  now  symmetric,  put  g,(^)  =  ff(-)' 
and  (2.1)  becomes  simply 


(2.3) 
with 

f2.3a' 


9[x)  +  9{y)  =  0  for  [x,  y)eSnD-,  X,  y  #  1 


g{z)  analytic  fox  \z\  <  1,  continuous  for  |;|  <  1,    z  ?^  1,  and 
lim(;-  l)g{z)  =  l-r4l,l)  =  1-p 


To  solve  this  consider  the  section  of  S  given  by  £  =  5  n  {y  =  z},  whose  equation  is 
\x\--r{x,x)  =  0.  Writing!  =  ^  +  tr;  we  find  that  £•  is  the  ellipse  {l+p)'{i~  ^)'  +  (t^W  = 
1  shown  in  figure  3. 


(.3,3 


Since  £  lies  in  the  unit  disk  except  at  a;  =  1,  (2.3)  implies  that 


(2.4) 


g{x)  +  g[x)  =  0  on  £",  x  ^  1  and 
lim  {x  —  l)g{x)  —  1  — p 


This  can  be  solved  using  conformal  mapping.   Let  <?i  be  the  conformal  mapping  from  the 
unit  disk  to  the  region  bounded  by  £,  normalized  by  ^(0)  =  -^— ■,  0(1)  =  1,  and  define 

h.{w)  =  g{4>[w)) 

_  (2.4)  is  now  easily  written  in  terms  of /i.  Let  P  denote  the  unit  circle,  D  the  open  and 
D  the  closed  unit  disk.  Note  that,  by  symmetry,  </)(ii')  =  <^(ti»)  and,  since  f  is  an  analytic 
curve,  4>  can  be  continued  analytically  to  |u)|  <  1  +  e  for  some  positive  e.  These  facts  imply 
that 


(2.5) 
(2.5a) 

(2.5b) 


h[w)  +  h{w)  =  0  on  r  \  {1},  with 
h  analytic  on  D,  continuous  on  £>  \  {1},  and 

lim(u,-l)/.(u,)  =  i^ 


Furthermore,  since  w  —  \/vj  on  F,  (2.5)  can  be  replaced  by 

(2.6)  h[w)  ^h{\/w)  =  Q  our  \{\) 

A  simple  variant  of  Schwartz  reflection  now  allows  us  to  extend  h  to  a  meromorphic  function 
in  the  extended  complex  plane  C"   =  C  U  {00},  by  defining  h{w)  =  -h{l/w)  for  \w\  >  1. 


Then  h  is  analytic  in  both  the  interior  and  exterior  of  the  unit  disk,  and,  by  (2.6),  h  is 
continuous  across  T  \  {l}.  Morera's  theorem  ensures  that  h  is  analytic  in  C"  \  {l},  and 
(2.5b)  implies  that  the  isolated  singularity  at  1  is  a  simple  pole  with  residue  (1  —p)/<f>'(l). 
From  these  conditions  we  see  that  h(w)  =  ^^  +  D  where  C  =  (1  -  p)/<?!>'(l).  Putting  this 
in  (2.6)  gives  D  =  C/2,  so  we  get 

Let  Tp  =  <f>~^  be  the  conformal  mapping  from  the  interior  of  £  to  D.  Then 
(2.8)  g(x)  =  MV'(^))^^(l-p)V''(l)^j^j^) 

The  function  V'  is  given  explicitly  by  (Nehari  [6]) 


(2.9)  V(a')  =  V^sn  (  —  sin-M  ^-^x- 1  >  "■  >      .,-(!_     II 


TT  \    p  )      )  \p      y  p- 


where  sn  is  the  Jacobi  elliptic  function  with  nome  q,   K   =   K{q)  and  k   —   k(q).     Using 
V'(l)  =  1  and  (sn'(;))^  =  (l  -  sn-(r))  (l  -yt=sn-(2)),  we  find  that 


(2.10)  ^'(i).^(l_,)yi±^ 

Putting  (2.9)  and  (2.10)  into  (2.8)  thus  gives  the  explicit  solution  for  g.  From  (2.2) 

(2.11)  p(x,y)=^^~^^^/^~^^g(x)+9(y)) 

xy-  r(x,y} 


(2.12)  =  2v/r 


p 


K{l-k)         (x-l)(3/-l)         V-(x)V.(y)-l 
V  (V'(x)  -  l)(V'(y)  -  1)    xy-T{x,y) 


3.    Solution  of  the  functional  equation:   Automorphy. 

In  this  section  we  show  how  a  simplified  version  of  the  method  of  Flatto  &  McKean  [4]  can 
be  used  to  construct  the  solution  g  to  (2.3).  The  simplification  involves  i)  ignoring  for  the 
moment  questions  of  uniqueness,  and  ii)  assuming  that  g  is  in  fact  a  meromorphic  function. 
We  know  that  g  has  a  simple  pole  at  1  with  known  residue.  If  xy  —  r(x,  y)  =  0  then  by 
(2.3)  a  simple  pole  at  x  implies  a  simple  pole  at  y  and  determines  its  residue.  Starting  from 
the  pole  at  1  we  thus  get,  as  we  shall  see,  a  countable  set  of  poles.  The  calculation  of  their 
locations  and  residues  is  made  possible  by  introducing  uniformizing  coordinates  for  the  curve 
xy  —  7-(x,y)  =  0.  Having  found  the  poles,  we  construct  a  partial-fraction  expansion  having 
these  poles  and  residues  and  show  by  direct  computation  that  it  satisfies  (2.3).  Uniqueness 
can  be  seen  either  from  the  reformulation  (2.5)  of  the  last  section  or  by  a  direct  appeal  to 
the  uniqueness  of  the  equilibrium  solution  p,j  derived  from  g  [2]. 

First  note  that  since  xy  —  r[x,y)  is  real,  symmetric  and  hyperbolic,  we  can  write 

xy  —  r[x,  y)  =  6~{xy  —  1),  where 
X  =  ax  +  (3y  +  7,  y  =  ay  4-  /3x  +  7 

with  a,  P,  7,  6  real.  The  inverse  relations  are 

X  ~  ax  +  Py  +  7  y  =  ay  +  /3x  +  7,  with 

2\         1+p        /  2\         1+p         /  1+p 


Now  the  curve  S  becomes  simply  y  =  1/i.  Define  X(x)  =  x(x,  1/x),  Y{x)  =  y(x,  1/x)  = 
X{l/x),  and  G(x)  -  g{X{x)).  Then  putting  i  =  X{x),  y  -  Y{x),  (2.3)  becomes 

(3.1)  G(f)  +  G(l/5)  =  0 

Rather  than  worry  about  where  this  holds,  we  look  for  a  G  that  is  meromorphic  and  satisfies 

(3.1)  for  aU  X  9^  0. 

For  I  fixed,  there  are  (except  at  branch  points)  two  values  of  y  satisfying  xy  —  r(x,  y)  =  0. 
In  terms  of  the  parameterization  of  5  by  £,  these  points  are  given  by  x  and  l/(Ax),  where 

Q    _    1  +   x/l  -  p^ 

That  is, 

(3.2)  X{x)  =  X{l/{\x)) 

This  impUes  that  G{x)  =  G(1/(A5)),  and  from  (3.1)  we  now  get 

(3.3)  G(x)  +  G{\x)  =  0 

Note  that  x  —  y  —  1  corresponds  to  i  =  y  =  1,  so  G  has  a  simple  pole  at  i  =  1  with 
residue  (1-p)  (^U=i)  =  \/l  -  p".  We  assume  that  G  is  meromorphic  in  C\{0}.  Then 
(3.3)  implies  that  G  has  simple  poles  at  i  =  A'",  m  £  Z,  with  residues  (  — A)'"-^!  —  p^. 
Therefore  g  has  simple  poles  at 

6^  =  A-(A'")  =dA"-+/3A-"'+7 
with  residues 

dX  . 


a^  =  (-A)-yr^^(— u- 


=  (-l)"'v'^l-p'(aA'"  -/3A-"') 

Since,  by  (3.2),  6_^_i  =  6^,  the  distinct  poles  correspond  to  m  >  0.  We  note  for  later  use 
that  6o  =  1  <  6i  <  62  •  ■  • ,  with 

vp      /  p- 

The  simplest  candidate  for  g  with  these  poles  and  residues  is 


1        "^        / 
(3.4)  ^(x)=i     ^        -^ 


T7l=I  —  00 


flm        \  P  -    1 


X  -  bm         a  -h,ri  )  '  P+1 


The   1/2  compensates  for  each  pole  appearing  twice,   and  the  second  term  both  ensures 
convergence  and  forces  g{a)  =  0,  which  follows  from  setting  x  =  — 1  in  (3.1). 

Now  verifying  that  g[x)  +  g{y)  =  0  for  all  (z,j/)  G  5  is  equivalent  to  showing  that 
G(i)  +  G(l/i)  =  0  for  all  i  7^  0,  where  G{x)  =  g(X(x)),  and  this  is  now  just  a  matter  of 
algebra: 

^     =G(i) 


V^ 


P- 


^X^(-i)'"(5A--/5a- 


1 

+ 


Q(i  -  A"") +/3(x-i  -  A-">)       qA"- + /9A-'" -f  7  -  a 


.00  . 

-y"(-i)'"(A'"  -  A-™-M  ( ^- + 


(x  -  A"')(l  -  x-iA-"*-!        A™  +  A-"'-!  +  7/d  -  a/a 
8 


Replacing  x  by  l/x  and  m  by  — m  gives  a  simOar  expression  for  G(l/x).   Adding  these  we 
get 

^         (G(i)  +  G(l/i)) 


yr^ 


T 


where  the  second  series  is  independent  of  x.  Therefore  we  get 


0 


(G(x)  +  G(i/f ))  =  ^ E(-ir  1  _  .-i\-n>.:  -  ,  _  ,-u-^-i  +  -"-^ 


Taking  every  other  term  in  the  series  gives  a  pair  of  convergent  telescoping  series  whose  sum 
is  0.  Therefore  G{x)  +  G{l/x)  is  constant,  and  setting  i  =  — 1  we  see  that  the  constant  is 
0.  Thus  g(x)  —  g(x)  is  the  unique  solution  to  (2.3). 

We  have  shown,  with  some  work,  that  g{x)  is  given  by  (3.4).  However,  we  could  have 
got  the  partial  fraction  expansion  of  G{x)  with  much  less  effort.  We  do  this  now,  since  for 
determining  the  moments  of  {pij},  G{x)  is  actually  eaisier  to  work  with  than  g(x),  as  we 
shall  see  in  the  next  section. 

G(x)  must  satisfy  (3.3)  for  x  ^^^  0,  with  the  additional  conditions  that 

(i)   G  is  meromorphic  for  x  ^  0,  with  simple  poles  only  at  i  =  A'",  m  £  Z;  the  residue 
at  1  is  \/l  —  p'. 

(ii)    G(— 1)  =  0,  obtained  by  setting  x  =  — 1  in  (3.1). 
Furthermore,  it  is  easy  to  see  that  this  determines  G  uniquely:    if  Gj,  Go  satisfied  these 
conditions  then  Gi  —  G2  would  be  analytic  in  C  \  {0}.  Therefore  Gi  -  G2  is  bounded  in  the 
annulus  1  <  |x|  <  A,  which  by  (3.4)  implies  Gi  —  G2  bounded  in  C,  hence  constant;  and  by 
condition  (ii)  the  constant  is  0.  Now  it  is  easy  to  write  down  the  expression  for  G: 

(-A)-     ,    (-A)-^ 


(3-5)  G(x)  =  ^r^-  Y,  ihr^^ 


m—  —  00 


1  +  A'^ 


Verifying  (i)  and  (ii)  is  most  easily  done  by  working  with  G'(i);  it  is  immediate  that  G'(x)  + 
AG'(Ax)  =  0,  from  which  G(x)  +  G(Ax)  is  constant.  Clearly  G(-l)  =  0,  and  G(-A)  =  0 
follows  with  a  Little  algebra,  so  the  constant  is  0,  and  the  conditions  are  verified. 

4.   Exact  and  asymptotic  results. 

Consider  the  expression  for  p(i,  y)  in  (2.11).  ^(x)  +  g{y)  is  meromorphic  in  each  variable 
and  vanishes  on  the  algebraic  curve  xy  —  r{x,y)  =  0  when  |x|,|y|  <  1.  Therefore  by  the 
principle  of  permanence  it  vanishes  everywhere  on  xy  —  r(x,  y)  —  0,  and  thus  the  only  poles 
of  p(x,j/)  are  those  of  g{x)  +  g{y),  except  for  the  poles  at  i  =  1  and  y  =  1,  which  are 
cancelled  by  the  factor  (x  —  l){y  —  !)• 

The  importance  of  this  observation  is  that  the  locations  of  the  poles  ofp(x,y),  considered 

as  a  function  of  x,  are  independent  of  y,  and  vice  versa.   This  means  that  the  asymptotic 

behavior  of  the  coefficients  p,j  is  determined,  as  in  the  univariate  case,  by  the  locations  of 

the  smallest  poles.   Specifically,  nunihor  the  poles  so  that  1  <  |Ai  |  <  jAo  |  <  ■  ■  • .   (Tiu-sc  arc 

the  same  with  respect  to  both  x  and  y  by  symmetry.)   Then  we  can  expand  p  successively 

in  X  and  y  as 

Arjy)        A.jy) 
P(-T,  y)  =  ;;-  +  r-  +  >f-\,,\, 

X   -     A\  X   "    A-> 


where  Rxfj.  denotes  a  function  regular  for  |i|  <  |A|,  \y\  <  |/Li|,  and  then 

2 

My)  =  J2  —^  +  ^--^= 


gives 


from  which 


('^»)=i:g,,-.;;;;-.,)^^>.>.^^>.>. 


with 


p_  =   Y,  A,,xr-'^r-'  +  O  ((IA3I  -  e)-"|Ai|—  +  |A,r-(|A3|  -  e)  — ) 


7l,j  =    lim     lim  (i  -  A,)(y  -  A_,)p(x,t/) 

X— 'A,  y— *Aj 


We  readily  see  from  (2.11)  that  An  =  0.  Therefore  the  asymptotic  behavior  of  the  equilib- 
rium probabilities  is  given  by 


for  J,  J  — ■  00,  where  A  =  A12  —  A2i- 

From  the  partial-fraction  expansion  (3.4)  of  g  in  section  3  we  find  that 

Ai  =bi  =QA-H/3/A-f  7  =  (--  1 

Vp 

,       -      ,  /         4        4        16        16 

A2  =  62  =  aX^  +  /3/A-  +7=1+ ^-  —  +  — 

V        p      p-       p"^       p^ 

Using  (2.11)  we  find  that 

yl  =  (Ai  -  1)(A2  -  l)ResgU,  (-^{^y  -  '■(^.  y))l(A.,A,)  j 

=  ^(l-p)=(4-p^) 

It  is  interesting  to  compare  the  asymptotic  behavior  of  p,j  to  that  of  the  equilibrium 
distribution  for  the  total  length,  N  =  I  +  J  of  the  two  queues,  pj,°'  =  X^,  +  ,=nPij.  and  to 
the  marginal  distribution  for  the  length  of  a  single  queue,  p^'"^  =  X^,  P.j-  The  generating 
functions  for  these  are,  respectively,  p(i,  z)  and  p(i,  1).   (2.11)  gives 

(4-2)  Pi-,-)=,^      %~'\^,9^^^ 

(1 -p-)(x-f  j^) 

from  which  we  find  (noting  that  the  pole  at  —  JT^  is  cancelled  by  a  zero  of  g)  that 


,„.       4(l-p)(2+p) 
Pn     ~(2_p)(2-p=)    1 


-A,- 

2-p'] 
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while  from  (1.3)  we  find 

(4.3)  p(z,l)=-i— - 

Ai    —  X 

SO  that 

(4.4)  p=""«  =  ^f^^^r 

The  non-independence  of  I  and  J  is  reflected  in  the  fact  that  p^'"^  and  pjj°'  ate  both 
asymptotically  proportional  to  A^",  whereas  if  they  were  independent  then  p5,°'  would  be 
asymptotically  proportional  to  nAJ"".  This  in  turn  reflects  the  fact  that  it  is  relatively 
unlikely  for  both  queues  to  be  long,  as  can  be  seen  by  comparing  p,j  (4.1)  with  p^'"^  '  Pj'"* 
(4.4).  For  instance  we  see  that  p.j  goes  like  (AiA2)~',  while  (p, )^  goes  like  (Aj)~'.  This  will 
also  be  seen  when  we  compute  covar(7,  J)  below. 

Now  consider  the  probability  poo  =  p(0,0)  that  both  queues  are  empty.  This  is  most 
easily  evaluated  using  the  expression  (2.12)  for  p[x,y)  in  terms  of  the  conformal  mapping 
function  t/'.  and  we  find  from  (2.9)  that 

Poo=— (1+P)^^  =  (1-P)-^^  =  (1-V^)  = 

TT 

where  K,  k  depend  on  p  via  the  nome  q.  For  instance,  when  p  =  1/2  we  get  (using 
Abramowitz  &;  Stegun  [l]  Table  17.3  to  obtain  m  =  k'  and  then  Table  17.1  for  K)  that 
Poo  =  .780.  As  p  — »  0,  Poo  — »  1  and  we  find  that  to  leading  order  (using  [l,  17.3.21,  17.3.22]) 

POO  =  1  -  ^p-  +  O(p^) 

Let  e(p)  be  the  probability  that  at  least  one  queue  is  empty.  Then  e(p)  =  2p(l,  0)  —  poo. 
For  example,  e(l/2)  =  .998.  From  (4.3),  p(l,0)  =  (1  -  p)/(l  -  p/2)-  .  Expanding  poo  to  6th 
order  in  p  via  the  formulas  cited  in  [1]  gives 

e(p)  =  1  -  ^P'  +  O(p') 

We  now  consider  the  problem  of  using  the  generating  function  p{x,y)  to  obtain  moments 
of  the  the  distribution  {pij}.  First  recall  that  because  of  the  form  of  (1.2)  we  had  the 
marginal  generating  function  p*'"8(x)  =  p(x,  1)  from  the  start  ((1.3)  and  (4.3)).  Therefore 
we  have  the  moments  £(/'^),  which  are  by  symmetry  =  E{J^).  In  particular, 

E[I)-       ^' 


var(7) 


4(1  -p) 
P^(p-2)^ 
16(1 -p)  = 


Thus  the  first  nontrivial  moment  which  p[x,y)  can  provide  is  E[JJ)  =  ^^^  p(l,  1),  which 
also  provides  the  variance  of  the  total  queue  length  N  =  I  +  J ,  since 

2 

(4.5)  E{N)  =  2E{I)  =        ^ 


2(1 -p) 
4.6)  var(Af)  =  2var(/)  +  2covar(7,  J) 
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In  fact,  it  is  easiest  to  work  with  the  generating  function  for  N,  p^°\x)  =  p(x,x),  given 
by  (4.2),  and  use  ^p'°'(l)  =  E[N{N  -  1)].   First  write 

o  1  p  —  1 

(4.7)  p'°'(z)  =  — ^ ;i(i),      where  a  =^^ -,      &nd  h{x)  =  {x  -  l)g{x). 

1  —  p'  X  —  a  p  +  i 

Note  that  h  is  analytic  at  1  and  h(l)  =  1  -  p.  The  moments  of  A^  are  thus  determined  by 
the  derivatives  of /i  at  1,  which  could  in  principle  be  determined  from  (2.8)  or  (3.4)  for  g. 
However,  in  view  of  (3.5),  it  is  simpler  to  work  in  the  transformed  variable  z  =  X{x)  of 
section  3.  Write 

(4.8)  H{i)  =  h{X{x))  =  ^l^^;\x  -  l)g{X{i)) 

X  —  1 

where  by  (3.5), 

(-A)'"         (-A)'" 


—  oo 

Substituting  this  and  X{x)  =  ax  +  P/x  +7  into  (4.8)  gives 

/  J  (  —  \\^         (  —  A) 

(4.9)  H[i)  =  Vl  -  p=(Q  -  /3/5)     1  +  :^(i  -  1)  +  (5  -  1)  E  k: 


+ 


i  -  A"'        1  +  A" 


Differentiating  (4.8)  and  using  X(l)  =  1  gives 


(4.10)  h'(\)  =  {X'{\))-'H'{\) 

h"{\)  =  {x\\)y'{H"{\)-h'[i)x"{\)) 

Using  (4.9)  we  find  that  h'(\)  =  1/2,  and 

if"(l)  =  -yr^^[/3+2(d-/3)  V      ^~^^' 


1  -x^y-  I 

which  gives 

'^"(^)  =  2(r^-^(^  +  ^)S(r^ 

Now  from  (4.7)  we  have 

which  yields 

dx^P      ^'^  2(1 -p)=  ^l-p  A.^(1-A™)2 

Smce  ^p"^'(l)  =  E{N^)  -  i;(A^),  we  get  (usmg  (4.5)) 


H.H)  var(Ar)  =  ZP!g±^P^_,i±^y- 


i-^Y 


Ail-pY-  i-p^Z-^Ji-A-n)^- 
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and  by  (4.6) 


(4.12)  covar(J,  J)  = 


16(1 -p)2         i_pZ-.^(i_A-)2 


Although  this  is,  in  the  words  of  Flatto  and  McKean,  "quite  unrevealing",  we  can  use 
it  to  find  the  asymptotic  behavior  of  covar(/,  J)  for  p  — •  1  ("heavy  loading").  First  note 
that  V^^o  (|:^L"j.  =  2  2^(-l)-A-"'(l  -  A-*")--.  Second,  A-i  ~  1  -  2V2^T^  for 
p  — *  1.  Furthermore,  since  ^^(-1)'"A"'"(1  -  A~'")~^  is  an  alternating  series  with  terms 
decreasing  in  absolute  value  to  0  we  can  take  limits  term  by  term  to  obtain 

L:=Um(l-p)f  <-^)"'^""'=f  izir 


1 


Now  y]-  ^lp  ^  2y:r  ^  -  Er  ^^  =  -(1/2)C(2),  wh.ch  gives  i  =  -:rV96  and 

Y^      (-A)"*  2L     _       -TT^ 

A.^(l_A-)2  ~  i_p-48(l-p)'      P^^ 

Substituting  in  (4.12)  gives 

27r-  -  21  -0.0263 


(4.13)  covar(/,J) 


48(1 -p)2        (l-p)=' 


Thus  the  two  queue  lengths  are  negatively  correlated,  at  least  for  p  —  1.  This  makes 
sense  intuitively  since,  at  each  time  step,  one  queue  can  lengthen  only  if  the  other  shortens 
(see  fig.  2),  so  that  if  7  is  large,  J  is  likely  to  be  small. 
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